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ABSTRACT 

We present the results of a program to monitor the four-image gravitational lens 
B1608+656 with the VLA. The system was observed over a seven month period from 
1996 October to 1997 May. The 64 epochs of observation have an average spacing of 3.6 d. 
The light curves of the four images of the background source show that the flux density of 
the background source has varied at the ~5% level. We measure time delays in the system 
based on common features that are seen in all four light curves. The three independent 
time delays in the system are found to be At B A = 31 ±7 d, Atsc — 36 ±7 d, and AtBD = 
76li d at 95% confidence. The uncertainties on the time delays are determined by Monte 
Carlo simulations which use fake light curves that have the characteristics of the observed 
light curves. This is the first gravitational lens system for which three independent time 
delays have been measured. A companion paper presents a mass model for the lensing 
galaxy which correctly reproduces the observed image positions, flux density ratios, and 
time delay ratios. The last condition is crucial for determining Ho with a four-image lens. 
We combine the time delays with the model to obtain a value for the Hubble constant 
of H = 59±7 km s" 1 Mpc" 1 at 95% confidence (statistical) for (Q m Ma) = (1,0). 
In addition, there is an estimated systematic uncertainty of ±15 km s -1 Mpc -1 from 
uncertainties in modeling the radial mass profiles of the lensing galaxies. The value 
of H Q presented in this paper is comparable to recent measurements of H from the 
gravitational lenses 0957+561, PG 1115+080, B0218+357, and PKS 1830-211. 

Subject headings: distance scale — galaxies: individual (B1608+656) — gravitational 
lensing 
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1. Introduction 

Even before the discovery of the first gravitational 
lens system, a technique for using gravitational lenses 
to measure the distance scale of the universe had been 
developed (Refsdal 1964). The technique requires a 
lens system in which multiple images of the back- 
ground source are formed. A "map" of the geodesies 
along which the light travels to form the images is con- 
structed and used to predict the differences in light 
travel times along the geodesies. If the background 
source is variable, these time delays can be measured 
as each image varies in turn. The ratios between the 
observed and predicted delays give the Hubble con- 
stant in the assumed world model (Om, ^a). The use 
of gravitational lenses for determining H has ma- 
jor advantages over traditional "distance ladder" ap- 
proaches. First, the technique gives a direct estimate 
of Ho at cosmological distances, where the effects of 
peculiar velocities are minimal. Second, this mea- 
surement of H is obtained in one step, without the 
propagation of errors inherent in the distance ladder 
approach. 

With the discovery that the "double quasar" 0957+561 
A,B was a lens system (Walsh ct al. 1979), the ef- 
fort to use lenses to measure H began in earnest. 
For many years the effort was hindered by both the 
paucity of known lens systems and the difficulty in 
measuring time delays in the systems. In fact, an 
unambiguous time delay has only recently been mea- 
sured in 0957+561 in spite of over ten years of in- 
tensive monitoring (Kundic et al. 1995; Kundic et 
al. 1997; Oscoz ct al. 1997; Haarsma et al. 1999). 
Another of the earliest known lenses, PG 1115+080 
(Weymann et al. 1980), has also just produced mea- 
surable time delays (Schechter et al. 1997). However, 
we may have entered a new era for time delay mea- 
surements from gravitational lenses. One reason for 
this is the accelerated rate of discovery of new lenses 
from systematic radio surveys. The Jodrell-VLA As- 
trometric Survey (JVAS; Patnaik et al. 1992; Browne 
et al. 1998; Wilkinson et al. 1998) has produced 6 
new lenses, and time delays have been measured for 
one of them (B0218357; Biggs et al. 1999). The on- 
going Cosmic Lens All-Sky Survey (CLASS; Myers et 
al. 1999) has found 12 new lenses since it began in 
1994. This paper reports the first measurement of 
time delays from a CLASS lens. 

The CLASS project is a large search for gravita- 
tional lenses with the VLA, with an explicit goal of 



finding lens systems which can be used to measure 
H Q . The gravitational lens B1608+656 (RA: 16 09 
13.956, Dec: +65 32 28.971, J2000) was observed in 
the first phase of CLASS and was immediately rec- 
ognized as a lens system. The radio discovery image 
shows four unresolved components in a typical lens 
geometry (Myers et al. 1995; see Fig. 1 for a map 
of the system). The system was also discovered in a 
search for gigahertz-peaked spectrum sources and was 
found to be the lensed core of a classical radio double 
source (Snellen et al. 1995). Further investigations of 
the system have provided data crucial for using the 
system for measuring Hq . We have measured the lens 
redshift (zg = 0.630; Myers et al. 1995) and source 
redshift (z s = 1.394; Fassnacht ct al. 1996). Optical 
and infrared observations taken with the Hubble Space 
Telescope reveal that the background source is being 
lensed by a pair of possibly merging galaxies (Jackson 
et al. 1997). The positions of the lensing galaxies rel- 
ative to the lensed images are important constraints 
on models of the lensing potential (Koopmans & Fass- 
nacht 1999; hereafter Paper 2). The HST images also 
show arcs due to the lensing of stellar emission from 
the background source. These arcs could be used as 
further constraints of the lens model. 

Flat-spectrum cores of radio galaxies such as the 
lensed object in B1608+656 are often variable. To 
test for variability in the B1608+656 background 
source, we made several observations of the system 
with the VLA, separated by time-scales of months. 
These data showed that the flux density of the back- 
ground source varied by up to 15%. The variability 
makes B1608+656 an excellent candidate for a dedi- 
cated monitoring program to determine time delays. 
This paper presents the results of VLA monitoring 
from October 1996 to May 1997. These observa- 
tions have resulted in the measurement of the three 
independent time delays in the B 1608+656 system 
and a subsequent determination of H n . The Hubble 
constant is expressed as H a = 100 h kms~ 1 Mpc~ 1 . 
Throughout this paper we assume (I^m, ^a) = (1, 0). 
The effect of varying the cosmological model is treated 
in Paper 2. 

2. Observations 

We observed B1608+656 between 1996 October 10 
and 1997 May 09, during which time the VLA was in 
the A, BnA, and B configurations. The 64 epochs 
were separated, on average, by 3.6 d. The obser- 
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vations were carried out at 8.5 GHz, giving angular 
resolutions ranging from 0'.'25 to 077 in the different 
array configurations. The observations are summa- 
rized in Table 1. The typical observation is 60 min 
long and includes scans on B1608+656, a flux calibra- 
tor (3C286 or 3C48), a phase calibrator (1642+689) 
chosen from the VLA calibrator list (Perley & Taylor 
f999), and two secondary flux calibrators (f 634+627 
and 1633+741). The secondary flux calibrators are 
nearby steep-spectrum sources that arc not expected 
to vary over the time-scales of the observations. We 
observe these sources to determine corrections for er- 
rors in the absolute flux calibration from epoch to 
epoch. The basic observing pattern is: 

• 1642+689 (1 min on source) 

• 1634+627 (1 min on source) 

• B1608+656 (4-6 min on source) 

• 1633+741 (2 min on source) 

A typical 60 min observation begins with a 9 - 
10 min scan (including slew time) on the flux cal- 
ibrator, contains three repetitions of the basic pat- 
tern on B1608+656, and ends with scans on 1642+689 
(1 min) and the flux calibrator again (^3 min). For 
the few observations that are 30 min or 90 min in 
length, the number of repetitions of the basic pattern 
is altered. 

3. Data Reduction 
3.1. Calibration 

The data reduction is separated into two major 
steps, calibration and mapping. The data for each 
epoch are calibrated using standard routines in the 
NRAO data reduction package AIPS. Before calibra- 
tion, the data quality for all sources is assessed and 
bad points are flagged with the EDITA and TVFLG 
tasks. Both of the flux calibrators are heavily resolved 
by the VLA in A configuration at 8.5 GHz. Hence, 
we cannot treat the calibrators as point sources with- 
out limiting the number of baselines that can be 
used to calculate phase and gain solutions. In or- 
der to increase the number of baselines available for 
the calculations, we create models of 3C286 and 
3C48 that incorporate the extended emission from 
the sources. We combine observations from several 
epochs with the DBCON task in AIPS. The resulting 



data sets have excellent (u, v)-plane coverage, from 
which we can make high dynamic-range maps. The 
mapping, which is performed in the DIFMAP package 
(Shepherd 1997), consists of alternating iterations of 
CLEANing (Hogbom 1974) and self-calibration. The 
final lists of CLEAN components are read back into 
AIPS and serve as the calibrator models. 

The procedures for the phase calibrator are simpler 
because the emission from 1642+689 is dominated by 
an unresolved component. Thus, the assumption that 
1642+689 is a point source leads to adequate phase 
and gain solutions. These calibration solutions are 
applied to the B1608+656, 1633+741, and 1634+627 
data. 

3.2. Source Maps 

We map the data and determine flux densities us- 
ing the DIFMAP package. We do not expect to see 
any structural changes over the course of the obser- 
vations; only changes in flux densities should be ob- 
served. To treat the data from each epoch in a uni- 
form fashion and to shorten the mapping procedure, 
we create models of the observed source structures 
from high dynamic-range maps. For epochs with 
noisy data, these models are needed to fix the loca- 
tions of the regions of low surface brightness emission, 
which otherwise would not be well-constrained by the 
data. 

To make the high dynamic-range maps of each 
source, we combine 13 high-quality data sets from A 
and B configuration observations with the DBCON 
task. The combined data sets, which have excellent 
(u, v) coverage, are then mapped in DIFMAP. All 
the maps are made with natural weighting. The sec- 
ondary flux calibrators, which have significant emis- 
sion from extended structures, are mapped by using 
an iterative cycle of CLEANing and self-calibration. 
Both phase and amplitude self-calibration are used. 
The models for these sources consist of the final 
lists of CLEAN components. The emission from 
B1608+656 is dominated by the four unresolved im- 
ages of the background source (Fig. 1). Hence, in- 
stead of CLEANing the data, we assign point-source 
model components to the four images of the back- 
ground source. We then use the DIFMAP modelf it 
function, which varies the component positions and 
flux densities to obtain the best fit to the (w, v)-plane 
visibilities. The model fitting iterations are alter- 
nated with phase and amplitude self-calibration. In 
later rounds of the model fitting, several nearby weak 
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Table 1 
Observations 







Array 


otan _l lmc 


ttot 


*1608 




JipOCll 


1V1 J LJ — O U U U U 


Configuration 


(^o l ) 


(min) 


(min) 


Comments 


1996 Oct 10 


366 


D^A 


16:30 


60 


27 




1996 Oct 12 


368 


D^A 


15:30 


30 


5 




1996 Oct 16 


372 


D^A 


21:00 


60 


21 




1996 Oct 18 


374 


A 


17:30 


60 


21 




1996 Oct 20 


376 


A 


01:30 


60 


21 


Elcv. < 30°. Thunderstorms 


1996 Oct 23 


379 


A 


17:30 


30 


6 




1996 Oct 26 


382 


A 


14:30 


60 


18 


T sys > 100K . 


1996 Oct 29 


385 


A 


22:30 


30 


5 


Elcv. < 30°. Ousting winds 


1996 Oct 31 


387 


A 


20:00 


60 


18 




1996 Nov 01 


388 


A 


18:30 


60 


18 




1996 Nov 03 


390 


A 


18:00 


60 


21 


Wind > 10 m/s 


1996 Nov 07 


394 


A 


07:00 


30 


6 


Elev. < 30° 


1996 Nov 08 


395 


A 


13:00 


60 


21 




1996 Nov 11 


398 


A 


13:00 


60 


18 




1996 Nov 12 


399 


A 


07:30 


30 


8 


Elcv. < 30° 


1996 Nov 15 


402 


A 


17:00 


60 


21 


High wind gusts 


1996 Nov 18 


405 


A 


15:30 


60 


21 




1996 Nov 23 


410 


A 


17:00 


60 


20 


Wind > 10 m/s 


1996 Nov 27 


414 


A 


21:00 


60 


20 


1996 Dec 01 


418 


A 


14:00 


60 


20 




1996 Dec 05 


422 


A 


14:30 


30 


7 




1996 Dec 07 


424 


A 


15:00 


30 


5 




1996 Dec 10 


427 


A 


17:00 


60 


20 


Wind > 10 m/s 


1996 Dec 15 


432 


A 


17:00 


90 


34 


1996 Dec 20 


437 


A 


19:15 


30 


5 




1996 Dec 23 


440 


A 


17:30 


90 


34 




1996 Dec 24 


441 


A 


15:30 


30 


10 




1996 Dec 29 


446 


A 


15:00 


60 


20 




1997 Jan 03 


451 


A 


16:30 


60 


20 


Wind > 10 m/s 


1997 Jan 07 


455 


A 


18:30 


90 


35 


1997 Jan 11 


459 


A 


16:00 


60 


20 




1997 Jan 16 


464 


A^BnA 


16:00 


60 


20 




1997 Jan 17 


465 


A^BnA 


09:30 


30 


8 


Elcv. < 30°. Flurries 


1997 Jan 20 


468 


A^BnA 


15:30 


60 


20 




1997 Jan 26 


474 


BnA 


14:00 


60 


20 


Rain 


1997 Jan 30 


478 


BnA 


13:00 


60 


19 




1997 Feb 02 


481 


BnA 


20:00 


60 


20 




1997 Feb 08 


487 


BnA 


19:00 


60 


19 




1997 Feb 13 


492 


BnA^B 


15:30 


60 


18 




1997 Feb 18 


497 


B 


18:00 


60 


19 




1997 Feb 23 


502 


B 


18:00 


60 


19 




1997 Feb 28 


507 


B 


18:00 


60 


19 


Snow storms. 


1997 Mar 03 


510 


B 


18:00 


60 


19 




1997 Mar 08 


515 


B 


16:00 


60 


18 




1997 Mar 14 


521 


B 


18:15 


75 


30 




1997 Mar 16 


523 


B 


17:45 


45 


12 




1997 Mar 21 


528 


B 


18:00 


60 


20 




1997 Mar 25 


532 


B 


21:30 


60 


19 


Wind > 10 m/s. Snowing. 


1997 Mar 30 


537 


B 


19:00 


90 


30 


1997 Apr 01 


539 


B 


17:30 


60 


20 




1997 Apr 05 


543 


B 


18:00 


60 


20 




1997 Apr 07 


545 


B 


18:00 


60 


20 




1997 Apr 11 


549 


B 


18:00 


60 


20 




1997 Apr 15 


553 


B 


18:00 


60 


20 




1997 Apr 19 


557 


B 


15:00 


60 


18 




1997 Apr 22 


560 


B 


18:00 


60 


20 




1997 Apr 26 


564 


B 


18:00 


60 


20 




1997 May 03 


571 


B 


18:00 


60 


20 




1997 May 08 


576 


B 


18:00 


90 


30 




1997 May 13 


581 


B 


18:00 


60 


20 




1997 May 17 


585 


B 


16:30 


60 


19 




1997 May 21 


589 


B 


18:00 


60 


20 




1997 May 23 


591 


B 


18:00 


60 


20 




1997 May 26 


594 


B 


20:00 


60 


19 
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sources are seen in the residual maps. These sources 
are included in the model for the last few iterations 
of the model fitting. The nearby sources can be seen 
in Fig. 2; their locations and flux densities are listed 
in Table 2. 

The first step in the mapping procedure at each 
epoch is to read in the data and perform a phase- 
only self-calibration against the model of the source. 
This procedure aligns the phase center of the obser- 
vation with that of the model and eliminates the need 
for many early steps of cleaning and self-calibration. 
After this point, the procedures used for B1608+656 
differ from those used for the secondary flux calibra- 
tors, as discussed below. 

3.2.1. B1608+656 

For each epoch, we determine the flux densities of 
the four lcnsed images in the B1608+656 system as 
follows. After the initial phase self-calibration step, 
the model described in the previous section is varied 
to find the best fit to the (u, w)-plane visibilities for 
that epoch. The component positions are held fixed 
and only the flux densities are allowed to vary. After 
several iterations of model fitting, another phase self- 
calibration is performed against the new model and 
more iterations of the model fitting are performed. 
At this point, the component flux densities and the 
RMS noise in the residual map are recorded in a log 
file (the "phasecal" data set). We then perform an 
overall gain calibration on the data, getting one gain 
correction per antenna for the observation. Typical 
gain corrections are on the order of 1 - 2%. Another 
round of model fitting is performed and the final com- 
ponent flux densities and RMS noise in the residual 
map are recorded in the log file (the "gscale" data 
set). We record the two separate data sets in case 
the gain calibration introduces any errors which may 
bias the subsequent analysis. All subsequent analy- 
sis is performed on both data sets and no significant 
differences are seen in the results. 

3.2.2. Secondary Flux Calibrators 




1 I I I I I I I I I I I I I I I I I I I I I c 

2 1 0-1-2 

Relative R.A. (arcsec) 

Fig. 1. — Map of B1608+656 from observation on 
1996 November 18. The contours are (—3, 3, 6, 12, 
24, 48, 96, 192, 384, 768) times the RMS noise level of 
0.035 mJy/beam. The map is made by fitting point 
source components to the (u, v) data and restoring 
with a 0'.'33xO'.'23 restoring beam. 



The steep-spectrum secondary flux-density calibra- 
tors contain significant extended emission and should 
not vary over the course of the observations. We ex- 
pect that any observed variations in flux density are 
due to errors in the absolute flux calibration, and as 
such can be expressed as an overall scaling of the 
model CLEAN component flux densities. That is, 
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Table 2 

Object Positions and Flux Densities 



Object 


a cr 


S (") a 


Ss.4 (mjy) b 


B1608+656 A 


+0.018 


+ 1.001 


34.3 


B1608+656 B 


+0.757 


-0.961 


16.8 


B1608+656 C 


+0.763 


+0.548 


17.4 


B1608+656 D 


-1.111 


-0.256 


5.9 


B1608+656 N 


+ 15.0 


+23.6 


0.6 


B1608+656 S 


-13.9 


-13.0 


0.5 


Source 1 


-34.7 


+ 167.0 


3.4 


Source 2 


-42.0 


-61.2 


0.8 


Source 3 


-178.0 


-43.5 


0.2 


Source 4 


-26.7 


+ 24.6 


0.2 



a Positions relative to map center at RA: 16:09:13.9530, 
Dec: +65:32:27.998 (J2000). 

Flux densities for B1608+656 components A - D arc 
the mean flux densities over the period of VLA monitoring 
(sec §3.4). 




200 100 -100 

Relative R.A. (arcsec) 



-200 



Fig. 2.— The field of B1608+656, showing the lens 
system and nearby radio sources. The object la- 
belled "1608" consists of the four components shown 
in Fig. 1. Objects "N" and "S" correspond to the 
northern and southern radio lobes of B1608+656 seen 
in low frequency maps (Snellen et al. 1995). 



the CLEAN component flux densities should not vary 
with respect to each other. The task of finding the 
overall flux density of these sources is thus simplified 
into finding the scale factor (fi) that, when multiplied 
by the component flux densities, gives the best fit to 
the data. To find the best-fit scaling for each data set, 
we create 11 scaled model files, based on the CLEAN- 
component models described above. The 11 files have 
/i ranging from 0.9 to 1.1 in steps of 0.02. Each of the 
scaled models is compared to the (u, w)-plane visibili- 
ties and a reduced \ 2 goodness-of-fit value is returned. 
We then fit a parabola to the points in the reduced- 
X 2 curve and find the value of \i which corresponds to 
the minimum reduced x 2 ■ This scaling gives the total 
flux density of the source at that epoch. 

3.3. Light Curve Editing 

Moore & Hewitt (1997), in their analysis of the 
15 GHz light curves of the gravitational lens MG 0414+0534, 
developed objective criteria to flag questionable data. 
They deleted from their light curves all points associ- 
ated with observations with the following conditions: 
the telescope elevations were less than 30°, the wind 
speed was greater than 10 m/s, or there was precipi- 
tation. We have noted all epochs satisfying their cri- 
teria in the "Comments" column of Table 1. How- 
ever, we are able to include many of these points in 
our analysis because observations at 8.5 GHz are less 
sensitive to the observing conditions than are obser- 
vations taken at 15 GHz. As such, we have only ex- 
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eluded epochs for which the data are severely affected 
by the observing conditions. This assessment is made 
by examining the light curves of the secondary cali- 
brator sources. All epochs for which the flux densities 
of the calibrators deviate by more than 15% from the 
mean value are deleted. Only two days are deleted 
after the application of this criterion: 376 and 382 
(MJD-50000). Note that epoch 376 satisfies two of 
the flagging criteria defined by Moore & Hewitt. At 
epoch 382, the system temperatures for all of the tele- 
scopes were in the range 100 - 200 K, as compared 
to the 30 - 50 K system temperatures measured for 
all other epochs. These high system temperatures 
may have resulted from the fact the subrcflcctors of 
several of the antennas had frozen prior to the obser- 
vation and had just thawed. The signal-to-noise ratio 
of the maps for epoch 382 were so low that no useful 
information could be extracted from them. The final 
edited light curves contain 62 epochs, with an average 
spacing of 3.7 d. Edited versions of the secondary flux 
calibrator and B1608+656 light curves are shown in 
Figures 3 and 4. 

3.4. Secondary Flux Calibration and Final 
Light Curves 

The normalized light curves of the secondary flux 
calibrators are shown in Fig. 3. As expected for non- 
variable sources, the light curves are close to con- 
stant. The scatter about the mean is small, with 
the exception of a few outlying points. All of these 
outlier points are due to bad observing conditions. 
The five worst points, for example, occur on days 
390, 402, 427, 474, and 532, all of which experienced 
high winds, precipitation, or both (Table 1). Because 
of these outliers, the formal calculation of a overes- 
timates the width of the error distribution. Conse- 
quently, we have estimated a using the interquartile 
range (IQR), which is less affected by the outliers. 
For a Gaussian distribution, a = 0.741 * IQR. For 
the secondary flux calibrators, we obtain a = 0.7% 
and a = 1.1% for 1634+627 and 1633+741, respec- 
tively. Gaussian distributions with these values of a 
are reasonable fits to the data (see Fig. 3). 

The light curves of the two secondary flux calibra- 
tors track each other extremely well, suggesting that 
the observed variations are due to errors in the ab- 
solute flux calibration of the data rather than any 
intrinsic variability of the sources. Because of this, 
we can use the light curves of the secondary flux 
calibrators to remove residual flux calibration errors 




g 0.9 



400 450 500 550 600 

MJD - 50000 




0.9 1 1.1 0.9 1 1.1 



Fig. 3. — Top- Light curves of the secondary cali- 
brators 1633+741 (open circles) and 1634+627 (filled 
boxes). Each curve has been divided by its median 
value over the length of the observations. The dotted 
vertical lines represent changes of array configuration. 
See Table 1 for the array configuration at each epoch. 
Bottom - The distribution of normalized flux densi- 
ties about the median values for 1634+627 (left) and 
1633+741 (right). 
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in the B1608+656 component light curves. We cre- 
ate a calibration correction curve by first normaliz- 
ing the light curves of the secondary flux calibra- 
tors and then averaging the two normalized fluxes 
at each epoch. We then apply the secondary flux 
calibration by dividing the B1608+656 component 
light curves by the calibration correction curve. The 
correction is effective in reducing the scatter in the 
B1608+656 light curves, as can be seen by compar- 
ing Fig. 4 to Fig. 5. The corrected component light 
curves are listed in Table 3. An electronic version 
of the light curve table may also be obtained from 
http : / /www . nrao . edu/~ cf assnac/1608f lux .tab. 

The errors on the final flux densities in Table 3 are 
a combination of additive and multiplicative terms. 
The additive uncertainty is well-approximated by the 
RMS noise in the residual map at the end of the 
model-fitting. The multiplicative term is indicative 
of how well the model-fitting procedure is able to find 
the "true" total flux density of a source. We estimate 
this error by calculating the ratio of the flux densities 
of 1634+627 and 1633+741, which we expect not to 
vary with respect to each other. The flux density ratio 
is not affected by errors in the absolute flux calibra- 
tion, so any scatter in the ratio can be attributed to 
errors in the model fitting procedure. From the scat- 
ter in the flux density ratio, we estimate the fractional 
error in the component flux densities to be 0.9%; the 
product of the component flux densities and this value 
gives the multiplicative uncertainty. The final errors 
on the flux densities are a combination in quadra- 
ture of the additive and multiplicative terms. For 
the three brightest components (A, B, and C), the 
multiplicative terms dominate; for component D the 
multiplicative and additive terms are comparable. 

The flux density of a lensed image of the back- 
ground source is Hi S u , where S v is the true flux den- 
sity of the background source and fii is the magni- 
fication associated with the i th (unresolved) image. 
Similarly, the lensed images will vary by /ii dS v if the 
flux density of the background source changes by an 
amount dS u . Hence, any variations in the B1608+656 
component light curves that are due to variations in 
the background source must have the same fractional 
amplitudes in all four curves. For this reason, the 
component light curves presented in Figs. 4 and 5 are 
divided by their mean flux densities to allow direct 
comparison of the fractional variations in the curves. 
The mean values used to normalize the light curves 
are 34.29, 16.79, 17.41, and 5.884 mJy for components 
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Fig. 4. — Light curves for the four images of the back- 
ground source in B 1608+656 before the secondary flux 
calibration has been applied. The light curves have 
been normalized by their mean values such that equal 
fractional variations in the flux densities of the com- 
ponents will have the same heights in the curves. 
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Fig. 5. — Same as Fig. 4 but after the application of 
the secondary flux calibration. 
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Table 3 

Final Component Light Curves 







Flux Density (mjy) 




MJD-50000 


A 


B C 


D 



366 


34 


.831 


± 





.309 


16. 


106 


± 





149 


17. 


602 


± 


0. 


162 


5, 


.822 


± 





.070 


368 


33 


.251 


± 





.307 


16. 


233 


± 


0. 


166 


17. 


269 


± 


0. 


175 


5, 


.736 


± 





.097 


372 


34 


.063 


± 





.308 


16. 


277 


± 


0. 


151 


16. 


643 


± 


0. 


154 


5, 


701 


± 





.063 


374 


34 


.003 


± 





.294 


16. 


532 


± 





146 


17. 


349 


± 


0, 


153 


5, 


760 


± 





.060 


379 


33 


.589 


± 





.307 


16. 


904 


± 





165 


17. 


042 


± 


0, 


166 


5, 


768 


± 





.085 


385 


32 


.972 


± 





.310 


17. 


002 


± 





176 


16. 


863 


± 


0. 


175 


5, 


695 


± 





.101 


387 


34 


.217 


± 





.317 


16. 


883 


± 


0. 


161 


16. 


831 


± 


0. 


160 


5, 


997 


± 





.070 


388 


33 


.381 


± 





.299 


16. 


884 


± 





155 


17. 


002 


± 


0. 


156 


6 


000 


± 





.066 


390 


34 


.845 


± 





.283 


17. 


365 


± 


0. 


149 


17. 


429 


± 


0. 


150 


5, 


.881 


± 





.073 


394 


32 


.744 


± 





.307 


16. 


811 


± 





174 


16. 


898 


± 


0. 


175 


5, 


727 


± 





.101 


395 


33 


.060 


± 





.297 


17. 


276 


± 





159 


16. 


715 


± 


0, 


154 


5. 


.763 


± 





.064 


398 


32 


.795 


± 





.295 


16. 


705 


± 


0. 


154 


16. 


995 


± 


0. 


157 


5, 


.755 


± 





.065 


399 


33 


.422 


± 





.304 


16. 


830 


± 





162 


16. 


661 


± 


0, 


161 


5. 


874 


± 





.082 


402 


33 


.927 


± 





.271 


16. 


943 


± 





139 


16. 


840 


± 


0. 


138 


6 


036 


± 





.060 


405 


33 


.846 


± 





.296 


16. 


740 


± 


0. 


150 


17, 


107 


± 


0. 


153 


5. 


910 


± 





.062 


410 


34 


.527 


± 





.310 


16. 


966 


± 


0. 


156 


17. 


294 


± 


0, 


159 


5. 


813 


± 





.064 


414 


33 


.878 


± 





.308 


17. 


096 


± 


0. 


160 


17. 


447 


± 


0. 


163 


5. 


813 


± 





.067 


418 


34 


.996 


± 





.315 


17. 


054 


± 


0. 


159 


17, 


632 


± 


0. 


164 


5. 


.807 


± 





.071 


422 


35 


.218 


± 





.315 


16. 


645 


± 


0. 


157 


17, 


431 


± 


0. 


164 


5. 


.765 


± 





.077 


424 


35 


.464 


± 





.323 


16. 


168 


± 


0. 


158 


17. 


324 


± 


0. 


168 


5. 


850 


± 





.084 


427 


35 


.909 


± 





.302 


16. 


613 


± 


0. 


144 


17, 


505 


± 


0. 


151 


5. 


.574 


± 





.061 


432 


34 


.811 


± 





.317 


16. 


423 


± 


0. 


152 


17, 


280 


± 


0. 


160 


5. 


663 


± 





.059 


437 


34 


.475 


± 





.308 


16. 


234 


± 





158 


17, 


368 


± 


0, 


167 


5. 


.791 


± 





.088 


440 


34 


.875 


± 





.310 


16. 


236 


± 


0. 


147 


17. 


532 


± 


0. 


158 


5. 


.483 


± 





.057 


441 


34 


.831 


± 





.315 


16. 


286 


± 


0. 


155 


17. 


203 


± 


0. 


163 


5. 


.525 


± 





.074 


446 


33 


.832 


± 





.304 


16. 


358 


± 





151 


17, 


405 


± 


0, 


160 


5. 


868 


± 





.064 


451 


34 


.444 


± 





.298 


16. 


320 


± 


0. 


146 


17. 


802 


± 


0. 


158 


5. 


.765 


± 





.066 


455 


33 


.230 


± 





.301 


16. 


271 


± 


0. 


150 


17. 


042 


± 


0. 


157 


5. 


825 


± 





.062 


459 


33 


.463 


± 





.298 


16. 


667 


± 





152 


16, 


946 


± 


0. 


154 


5. 


964 


± 





.064 


464 


34 


.451 


± 





.310 


16. 


851 


± 


0. 


156 


17. 


364 


± 


0. 


160 


5. 


979 


± 





.068 


465 


34 


.010 


± 





.309 


16. 


720 


± 


0. 


164 


17, 


825 


± 





173 


6 


049 


± 





.090 


468 


33 


.989 


± 





.305 


16. 


624 


± 


0. 


152 


16, 


867 


± 


0. 


154 


5. 


970 


± 





.063 


474 


33 


.015 


± 





.322 


16. 


290 


± 





172 


16, 


755 


± 


0. 


176 


5. 


974 


± 





.095 


478 


32 


.845 


± 





.295 


16. 


223 


± 


0. 


149 


16. 


945 


± 


0. 


156 


5. 


845 


± 





.064 


481 


33 


.523 


± 





.306 


16. 


735 


± 


0. 


155 


16, 


709 


± 


0. 


155 


6. 


015 


± 





.063 


487 


33 


.643 


± 





.306 


16. 


325 


± 


0. 


151 


17, 


410 


± 


0. 


161 


5. 


960 


± 





.062 


492 


34 


.166 


± 





.308 


16. 


528 


± 





152 


16. 


923 


± 


0. 


155 


5. 


820 


± 





.062 


497 


33 


.754 


± 





.304 


17. 


080 


± 





156 


16. 


961 


± 


0, 


155 


5. 


725 


± 





.060 


502 


34 


.051 


± 





.307 


16. 


427 


± 


0. 


150 


17. 


323 


± 


0. 


158 


5. 


879 


± 





.060 


507 


34 


.040 


± 





.308 


16. 


761 


± 


0. 


156 


17. 


495 


± 


0. 


162 


5. 


808 


± 





.067 


510 


33 


.924 


± 





.307 


16. 


779 


± 





154 


17, 


754 


± 


0. 


163 


5. 


.779 


± 





.061 


515 


32 


.624 


± 





.292 


17. 


156 


± 


0. 


156 


17. 


050 


± 


0. 


155 


5. 


804 


± 





.060 


521 


33 


.901 


± 





.307 


16. 


774 


± 


0. 


154 


17, 


356 


± 


0. 


159 


5. 


.743 


± 





.060 


523 


34 


102 


± 





.310 


16. 


648 


± 





157 


17, 


337 


± 





162 


5. 


.792 


± 





.069 


528 


34 


.548 


± 





.314 


16. 


691 


± 


0. 


155 


17. 


386 


± 


0. 


161 


5. 


582 


± 





.061 


532 


34 


.210 


± 





.295 


16. 


394 


± 





145 


17, 


069 


± 


0. 


150 


5. 


777 


± 





.061 


537 


33 


.749 


± 





.307 


16. 


636 


± 


0. 


153 


17. 


181 


± 


0. 


158 


5. 


.792 


± 





.058 


539 


34 


.388 


± 





.309 


16. 


958 


± 


0. 


155 


16, 


927 


± 


0, 


155 


5. 


826 


± 





.061 


543 


34 


.049 


± 





.306 


16. 


921 


± 





155 


17, 


185 


± 


0. 


157 


5. 


830 


± 





.062 


545 


33 


.959 


± 





.308 


17. 


031 


± 


0. 


156 


17. 


422 


± 


0. 


160 


5. 


903 


± 





.060 


549 


34 


.338 


± 





.311 


16. 


267 


± 


0. 


150 


17, 


507 


± 


0. 


161 


5. 


.758 


± 





.061 


553 


33 


.728 


± 





.306 


16. 


519 


± 


0. 


152 


17. 


374 


± 


0. 


160 


5. 


820 


± 





.061 


557 


33 


.004 


± 





.295 


16. 


755 


± 





152 


17, 


459 


± 


0, 


158 


5. 


.707 


± 





.060 


560 


33 


.241 


± 





.300 


16. 


622 


± 


0. 


152 


17, 


630 


± 


0, 


161 


5. 


.833 


± 





.061 


564 


33 


.049 


± 





.300 


16. 


850 


± 


0. 


156 


16, 


666 


± 


0. 


154 


5 


.805 


± 





.062 


571 


34 


.091 


± 





.309 


16. 


535 


± 





152 


17, 


653 


± 


0, 


162 


5. 


957 


± 





.063 


576 


33 


.332 


± 





.302 


16. 


748 


± 


0. 


154 


17. 


093 


± 


0. 


157 


5. 


.847 


± 





.059 


581 


34 


.055 


± 





.308 


15. 


971 


± 


0. 


147 


18, 


138 


± 


0, 


166 


5. 


768 


± 





.060 


585 


34 


.151 


± 





.307 


16. 


127 


± 





148 


17, 


287 


± 


0. 


158 


5. 


.732 


± 





.060 


589 


33 


.723 


± 





.307 


15. 


952 


± 


0. 


148 


17. 


538 


± 


0. 


162 


6. 


.051 


± 





.064 


591 


34 


.341 


± 





.312 


16. 


068 


± 





150 


17, 


214 


± 


0, 


160 


5. 


962 


± 





.066 


594 


33 


.130 


± 





.302 


16. 


309 


± 


0. 


152 


17, 


018 


± 


0. 


158 


5. 


993 


± 





.065 
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A, B, C, and D, respectively. All four light curves 
show variations in flux density at the ~5% level, and 
there are common features that can be seen in each. 
Each light curve shows a rise of ^5% in flux density, 
followed by a plateau lasting ~20 d and then a drop 
of -4%. 

4. Determination of Time Delays 

Models of the B1608+656 system predict that, if 
the background source is variable, the lcnscd images 
will vary in the order B — > A — > C — > D (Myers 
et al. 1995; Paper 2). This is, in fact, the behav- 
ior seen in the B1608+656 light curves. In order to 
determine time delays between component B and the 
other three components, we have used three statisti- 
cal methods: (1) smoothing and x 2 minimization, (2) 
smoothing and cross-correlation, and (3) dispersion 
analysis (Pelt et al. 1994, 1996, 1998). The methods 
are discussed below. All of the analysis has been per- 
formed on both the "phasecal" and "gscale" data sets 
(§3.2.1). No significant differences in the results are 
seen, so we present only the results from the "gscale" 
set in the subsequent discussion. 

4.1. Methods Using Smoothing/Interpolation 

The first two methods of determining the time de- 
lays from the B1608+656 light curves require that 
the observed data be transferred onto a regular grid. 
Some previous determinations of time delays (e.g., 
Kundic et al. 1997; Biggs et al. 1999) have accom- 
plished this transfer by linear interpolation of the 
data. However, the flux density variations seen in the 
B1608+656 images are small compared to the noise 
in the curves; therefore, linear interpolation can am- 
plify noise spikes. In contrast, smoothing reduces the 
effects of noise compared to the true variations, if the 
variations have typical time-scales significantly longer 
than the sampling interval. We smooth and re-sample 
the data by calculating the weighted mean of points 
within a smoothing window that is moved from the 
beginning to the end of the observations in regular 
steps. The step size is set to one day. 

We smooth each light curve with several different 
functions to avoid biasing the results by our choice 
of weighting function or window size. For complete- 
ness, we also include the results obtained from in- 
terpolating the data by piecewise linear interpola- 
tion. We use the following smoothing schemes: (1) 
boxcar-weighted mean, (2) triangle-weighted mean, 



(3) Gaussian- weighted mean, (4) boxcar- weighted mean 
with a variable-width smoothing window, and (5) 
triangle-weighted mean with a variable-width win- 
dow. In the last two schemes, the width of the 
smoothing window is varied such that the same num- 
ber of points are always included in the window. We 
use 3, 5, and 7 point windows for these methods. For 
the fixed- width window schemes, we smooth with win- 
dow widths of 5, 10, and 15 d. The Gaussian smooth- 
ing scheme uses values of a equal to 3, 5, and 7 d. In 
addition to multiplying by the value of the smoothing 
function, the data points are also variance-weighted. 
Thus, the overall weighting on a point at t — ti which 
is being used to calculate a weighted mean at t = tk 
is: 

■^g(U - t k ) 
T,jl-^9{tj -t k )] 

where g(t) is the value of the smoothing function (e.g., 
g(t) = 1 for boxcar smoothing), Oi is the uncertainty 
on the flux density at t = U (calculated in §3.4), and 
the sum is over the points used to calculate the mean. 
The variance of the weighted mean y at t — tk in the 
smoothed curve then becomes 



2 

a y,k 



1 



(Ei{^9(ti-tk)}) 



.£[£&-*>!!, (2) 



and (Jy : k is taken as the uncertainty in the smoothed 
flux density at that step. 

4-1.1. x 2 Minimization 

For the x 2 minimization technique, we compare 
two light curves, one of which is designated the "con- 
trol" curve and the other of which is designated the 
"comparison" curve. Both curves are smoothed, and 
then the comparison curve is multiplied by a scale 
factor /i so that its mean flux density is comparable 
to that of the control curve. After this scaling, the 
comparison curve is shifted in time by an amount t 
with respect to the control curve. We form a grid of 
(p,, t) pairs and calculate the x 2 statistic at each grid 
point. The step size for r used in the grid is 1 d and 
the minimum and maximum shifts are set to ±114 d, 
i.e., half of the total length of the observations. The 
amplitude scale factors are set as percentages of the 
scale factor fi that equalizes the mean flux densities 
of the two curves being compared. The scale factors 
used in the grid range from 0.9/xo to 1.1/Uo, in steps 
of O.OOl^o- 
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One result of the smoothing and interpolation per- 
formed on the input light curves is that the points 
in the interpolated curves are no longer independent. 
When computing the reduced x 2 statistic we estimate 
an effective number of independent points, N e ff, in 
the overlapping region between the shifted and un- 
shifted curves. This quantity is set to the number 
of points that would have been present in the over- 
lapping region if the original light curves had been 
regularly sampled at the average spacing of 3.7 d. 
This simplifying assumption should not significantly 
affect the quantity in which we are interested, which 
is the location of the minimum on the \ 2 surface. In 
fact, the ^-minimization method finds delays that 
match closely those found from the completely non- 
interpolative dispersion method (§4.2), indicating the 
robustness of the delays determined from this method. 

We repeat the x 2 minimization three times, once 
for each independent pair of light curves. The compo- 
nent B light curve is always taken as the control curve. 
In the three repetitions of the \ 2 minimization, the 
comparison curves are A, C, and D, respectively. The 
minimum delay is found by fitting a parabola to the 
points at the minimum of the gridded x 2 curve. Ta- 
ble 4 summarizes the results. Typical goodness-of-fit 
curves for the three pairs of light curves are shown in 
Figure 6. 

4-1-2. Cross- correlation 

For the cross-correlation calculations the compo- 
nent B light curve is once again taken as the control 
curve. Before the cross-correlations are performed, 
both the comparison and control curves are smoothed 
and then divided by their mean values over the obser- 
vations. The normalization by the mean value puts 
the fractional variations in the light curves at the 
same level. The cross-correlation functions are com- 
puted in the time domain, with the value at each lag 
calculated in the standard fashion: 



N X 20 




: BA 




BC 




:: BD 



-100 100 -100 100 -100 100 
Delay Delay Delay 

Fig. 6. — Plots of reduced \ 2 vs - lag from comparison 
of the light curves of components B and A (left), B 
and C (middle), and B and D (right). Each light 
curve is smoothed with a boxcar of width 10 d before 
the x 2 curves are calculated. The minima (vertical 
dashed lines) are at lags of 29.3, 36.2, and 75.3 d, 
respectively. 



r.ik = 



s 2 

_ s jk 



SjS k 



(3) 



where s 2 fc is the covariancc defined by 



„2 _. 



N-l Si a^.+a\ ( x ij x j)(xik Xk) 



N Li al+aj h 



(4) 



and N is the number of overlapping points (e.g., Bev- 
ington 1969). Note that in Xij or ov, , the first sub- 
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Table 4 

Time Delays from x 2 Minimization and Cross- Correlation 



Smoothing 


Smoothing 


x 2 


Minimizat 


ion 


Cross-correlation 




Function 


"Width" a 


At B A 


At B 


c 


At BD 


AtBA 


At BC At 


BD 


Boxcar 


5 


30.97 


35. 


96 


74.27 


28 


36 


74 




10 


29.25 


36. 


24 


75.32 


29 


37 


75 




15 


31.84 


36. 


10 


75.58 


32 


36 


75 


Triangle 


5 


29.16 


36 


07 


74.20 


28 


36 


74 




10 


29.35 


35. 


63 


75.27 


29 


35 


75 




15 


30.77 


36. 


03 


75.97 


30 


36 


76 


Gaussian 


3 


30.65 


35. 


93 


75.93 


30 


36 


76 




5 


31.55 


36. 


17 


76.70 


31 


37 


76 




7 


31.52 


36 


.37 


77.06 


31 


37 


77 


Variable- width 


3 


30.13 


35. 


.79 


76.58 


30 


36 


75 


boxcar 


5 


29.92 


36. 


.18 


77.32 


30 


37 


78 




7 


29.44 


36. 


.38 


77.22 


30 


39 


77 


Variable- width 


3 


30.52 


35. 


.84 


76.48 


30 


36 


76 


triangle 


5 


30.74 


36 


22 


77.03 


30 


37 


77 




7 


29.93 


36. 


.25 


77.21 


30 


37 


77 


Interpolation 




29.02 


36. 


09 


74.82 


28 


36 


74 



a "Width" is width of window for boxcar and triangle smoothings, a for Gaussian smoothing, 
and number of points in window for variable-width smoothings. 



Table 5 

Flux Density Ratios from \ 2 Minimization 



Smoothing 


Smoothing 








Function 


"Width" 


Sa/Sb 


Scl Sb 


Sd/ Sb 


Boxcar 


5 


2.0387 


1.0368 


0.3506 




10 


2.0383 


1.0378 


0.3508 




15 


2.0411 


1.0379 


0.3508 


Triangle 


5 


2.0387 


1.0368 


0.3505 




10 


2.0400 


1.0376 


0.3507 




15 


2.0409 


1.0379 


0.3511 


Gaussian 


3 


2.0407 


1.0378 


0.3511 




5 


2.0411 


1.0379 


0.3511 




7 


2.0410 


1.0379 


0.3508 


Variable-width 


3 


2.0406 


1.0383 


0.3511 


boxcar 


5 


2.0415 


1.0384 


0.3512 




7 


2.0413 


1.0382 


0.3512 


Variable- width 


3 


2.0407 


1.0383 


0.3511 


triangle 


5 


2.0414 


1.0384 


0.3512 




7 


2.0393 


1.0382 


0.3512 


Interpolation 




2.0394 


1.0374 


0.3506 
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script refers to a particular observation, and the sec- 
ond subscript refers to the name of the variable under 
discussion. Thus the i th observation of variable 

Xj. The variance, s 2 , is given by 



N- 



-1 Ei ^( x >j x jY 



V Si 



(5) 



N 



where the means are weighted means, 



E* 



(6) 



The cross-correlation calculations are repeated for 
all of the smoothing functions described in §4.1. In all 
cases, clear peaks in the cross-correlation curves are 
seen. Typical curves are shown in Fig. 7. The lags at 
which the peak correlation coefficients occur are given 
in Table 4. The average displacement between the 
cross-correlation and \ 2 minimization lags for each 
pair of curves is less than one day. 

4.2. Dispersion Analysis 

The dispersion analysis methods presented by Pelt 
et al. (1994, 1996, 1998) do not involve any interpo- 
lation of the component light curves. These methods 
thus have the advantage of avoiding effects introduced 
by the interpolation and smoothing associated with 
the methods discussed in §4.1. The dispersion analy- 
sis begins with the construction of a composite curve 
Cfe(tfe) from two input light curves Ai and Bj. Curve 
Ai is not modified, while curve Bj is scaled by a flux 
density ratio fj, and shifted in time by a delay r, i.e. 



Ai, 
/iBj, 



if tk = ti 

if tk ^t-j+r 



(7) 



(Pelt et al. 1996). We have used a grid of (/U,r) 
pairs to construct the Ck curves. Aside from chang- 
ing the spacing on the delay axis to 0.5 d, the grid 
limits and spacings are the same as those used in the 
X 2 minimization analysis presented in §4.1.1. The 
internal dispersion in each curve is calculated, and 
the grid point associated with the minimum disper- 
sion is recorded. In our analysis of the B1608+656 
light curves we calculate dispersions using the non- 
parametric D\ and the one-parameter D 2 2 statistics, 
where the notation is taken from Pelt ct al. (1996; 
note that these statistics are called D\ and D 2 , re- 
spectively, in Pelt et al. 1998). The D 2 dispersion is 






-100 100 -100 100 -100 100 
Delay Delay Delay 

Fig. 7. — Plot of correlation coefficient vs. lag from 
comparison of the light curves of components B and 
A (left), B and C (middle), and B and D (right). 
Each light curve was smoothed with a Gaussian with 
<j = 5d before the cross-correlations were calculated. 
The maxima (vertical dashed lines) are at lags of 31, 
37, and 76 d. 
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calculated using only immediately adjacent points in 
the composite curve, with the caveat that a pair of 
points only contributes to the dispersion if the two 
points are from different input curves. The D\ 2 dis- 
persion is similar, but uses all pairs of points that 
lie within a time range 5 of each other. For a de- 
tailed description of these estimates, see Pelt et al. 
(1996). Table 6 gives the results for the D\ analysis 
and the D\ 2 analysis with several values of the 6 pa- 
rameter. The results do not depend strongly on the 
value of S and are consistent within the errors (§5.2) 
with the results from the x 2 minimization and cross- 
correlation analyses. The dispersion spectra plotted 
in Fig. 8 are cuts through the (fi, r) grid at a value of 
fi corresponding to the minimum dispersion. 

4.3. Time Delays 

There is some scatter in the results presented in 
Tables 4 and 6, which is not surprising considering 
the low levels of variation and sparse sampling of the 
light curves. However, the scatter is small compared 
with both the length of the time delays and the un- 
certainties in the delays that we find from the Monte 
Carlo simulations in §5.2. For each of the three meth- 
ods used to find the delays (x 2 minimization, cross- 
correlation, and dispersion analysis) , we take the me- 
dian values of the delays in Tables 4 and 6. The three 
median values for each delay are then averaged to 
obtain AtsA = 31 d, Atsc — 36 d, and AtsD = 
76 d. The median flux density ratios, computed in 
a similar manner, are found to be Sa/Sb = 2.0418, 
Sc/Sb = 1.0376, and S D /S B = 0.3512. We shift the 
light curves by the mean delays and normalize them 
using the mean flux density ratios to create a com- 
posite light curve of the background source (Fig. 9). 
A composite curve constructed from the smoothed 
and interpolated component light curves is shown in 
Fig. 10. 

5. Monte Carlo Simulations 

5.1. Significance of Light Curve Correlations 

The variations seen in the B1608+656 light curves 
are not large, either in a fractional or absolute sense, 
compared to what has been seen in other lens sys- 
tems. The fractional variations seen in B0218357, 
0957+561, and PG 1115+080 are all two to three 
times larger than those seen in B1608+656 (Biggs et 
al. 1999; Kundic et al. 1995,1997; Schechter et al. 
1997), and PKS 1830-211 shows 50% variations in 



flux density (although the time delay measurement is 
based on smaller variations; Lovell et al. 1998). A 
sceptic might argue that the correlations between the 
B1608+656 light curves are not significant and could 
be duplicated by any set of light curves containing 
random scatter about a constant value. 

In theory, the value of the correlation coefficient 
can be used to assess the significance of the correla- 
tion. With the coefficient, r calculated as described 
in §4.1.2, the probability of obtaining a value > r 
from two uncorrelated curves of Gaussian-distributed 
random variables is 



P c (r,A) = 2 f P r (p,u)dp 

J\r\ 



where 



P r {r, v) 



1 rhV+l)/2] 
V5F r(i//2) 



(1 



r 2j(v-2)/2 



(8) 



(9) 



and v — N — 2 (e.g., Bevington 1969). However, 
the smoothing and interpolation performed on the 
sparsely sampled B1608+656 light curves makes the 
interpretation of the significance of the value of r 
complicated. The difficulty lies in assessing the num- 
ber of independent points in the curves at each lag. 
Monte Carlo simulations show that the number of in- 
dependent points cannot be estimated simply as the 
width of the overlap region divided by the width of 
the smoothing window (i.e., the number of smooth- 
ing windows in the overlap region). This calculation 
underestimates the number of independent points in 
the region. Because it is difficult to determine the 
significance of the correlation analytically, we per- 
form Monte Carlo simulations to find the significance 
empirically In the simulations, we calculate correla- 
tions between light curves consisting of randomly dis- 
tributed data. To create random light curves with the 
same distribution of flux densities as seen in the data, 
we simply randomize the time series for the compo- 
nent light curves while preserving the flux densities at 
their measured values. By randomizing the times at 
which the flux densities are measured, we destroy any 
possible correlations between the curves. Each light 
curve is randomized independently to avoid correla- 
tions at zero lag which are associated with measure- 
ment errors. 

The simulations are conducted with 3000 sets of 
randomized curves. Each set of light curves is pro- 
cessed in the manner described in §4.1.2 and produces 
three sets of correlation curves (B-A, B-C, and B-D). 
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Table 6 
Dispersion Analysis Results 



Statistic 


<5 


At BA 


At BC 




Sa/Sb 


Sc/Sb 


Sd/Sb 


n 2 




30 5 




36 5 


1 0372 


74 5 


V .oovo 


n 2 


3.5 


32.0 


2.0429 


35.5 


1.0361 


74.5 


0.3501 


n 2 


4.5 


31.5 


2.0429 


36.5 


1.0372 


74.5 


0.3502 


D 2 


5.5 


31.0 


2.0429 


37.0 


1.0372 


76.5 


0.3509 


D 2 


6.5 


31.0 


2.0429 


34.5 


1.0361 


76.5 


0.3509 


D 2 
^4,2 


7.5 


31.5 


2.0429 


35.0 


1.0372 


76.5 


0.3509 


D 2 

^4,2 


8.5 


32.0 


2.0429 


35.0 


1.0372 


77.0 


0.3512 


D 2 

^4,2 


9.5 


31.5 


2.0449 


34.5 


1.0372 


77.5 


0.3512 


n 2 

^4,2 


10.5 


31.5 


2.0449 


35.0 


1.0372 


77.0 


0.3512 


D 2 

^4,2 


11.5 


32.0 


2.0449 


35.5 


1.0372 


77.5 


0.3512 


D 2 

^4,2 


12.5 


32.0 


2.0449 


36.5 


1.0372 


78.5 


0.3512 


D 2 

^4,2 


13.5 


32.0 


2.0449 


35.5 


1.0372 


78.5 


0.3512 


D 2 
^4,2 


14.5 


32.5 


2.0449 


35.0 


1.0372 


78.5 


0.3512 


D 2 

^4,2 


15.5 


33.0 


2.0449 


34.5 


1.0372 


77.5 


0.3512 


D 2 

^4,2 


16.5 


32.5 


2.0449 


35.0 


1.0372 


79.5 


0.3516 


D 2 
^4,2 


17.5 


32.5 


2.0449 


35.0 


1.0372 


79.0 


0.3512 


ZJ 2 

^4,2 


18.5 


32.0 


2.0429 


35.0 


1.0372 


79.5 


0.3512 


D 2 
^4,2 


19.5 


32.0 


2.0429 


35.0 


1.0372 


79.5 


0.3512 


D 2 

^4,2 


20.5 


32.0 


2.0429 


35.0 


1.0372 


80.0 


0.3516 



a 
w 

s 






-100 100 -100 100 -100 100 
Delay Delay Delay 

Fig. 8. — Plot of dispersion spectra from comparison 
of the light curves of components B and A (left), B 
and C (middle), and B and D (right). The spectra 
were calculated with the D\ 2 method with <5 = 5.5d. 
The minima (vertical dashed lines) are at lags of 31, 
37, and 77 d. 
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Fig. 9. — Composite light curve constructed by nor- 
malizing the component light curves and shifting by 
the lags given in §4.3. The light curves are repre- 
sented by squares (component A), diamonds (B), tri- 
angles (C), and circles (D). For clarity of presentation, 
error bars represent only the additive component con- 
tributed by the RMS noise in the maps (see §3.4). 
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All values of the correlation coefficient are recorded. 
The distributions of the cross-correlation values ob- 
tained from the 10 d boxcar smoothing scheme are 
shown in Fig. 11. The empirical probabilities of ob- 
taining at least the observed peak values, which are 
indicated by the vertical dashed lines in the figure, 
from uncorrelated curves are all low, with P(\r\ > 
\r A \) = 2.8 x 1(T 4 , P(\r\ > \r c \) = 4.3 x 1CT 3 , and 
P(\ r \ > l r r>l) < 1-5 x 10~ 6 . It is even more unlikely 
that three pairs of randomized curves could produce 
three such anomalously high cross-correlation peaks. 
There can thus be no significant doubt that the cor- 
relations we measure are real. 

5.2. Uncertainties in Time Delays and Flux 
Density Ratios 

The time delay measurement uncertainties con- 
tribute directly to the error budget for measuring Ho 
with a gravitational lens. In particular, the fractional 
uncertainties in the time delays contribute a matching 
fractional uncertainty in H 0l i.e., 

(vh) delay _ ^At ^ 

Hq At 

We estimate the uncertainties in the time delays by 
performing Monte Carlo simulations of the observa- 
tions. In the simulations we assign time delays be- 
tween the pairs of light curves and then see how well 
we can recover the input delays. We also use the 
simulations to estimate the uncertainties in the flux 
density ratios, which are necessary for modeling the 
lensing potential (Paper 2). 

We produce fake curves with the same characteris- 
tics as our real data by smoothing the composite light 
curve (Fig. 9) with a 10 d boxcar filter. This normal- 
ized and smoothed curve is the master light curve for 
the simulations, representing the assumed true be- 
havior of the background source. The offsets between 
the points of the composite curve and the master 
curve are distributed as a zero-mean Gaussian with 
c = 0.014. Note that this is & fractional value since all 
of the light curves have been normalized to create the 
master light curve. The Gaussian distribution is used 
to generate the random offsets for the simulations. 
The appropriate rescaling of the random offsets is 
achieved through Equation 11. We note that smooth- 
ing procedure used in constructing the master curve 
destroys information on possible short-timescale vari- 
ations. Thus, if such short-term variations do exist, 
their effects will be attributed to measurement and 



calibration errors in this method. However, given the 
low level of variability of the background source dur- 
ing the observations and the sparse sampling, it is dif- 
ficult to distinguish any possible short-term variations 
from measurement error. Thus, we accept o = 0.014 
as an indication of the measurement error, but with 
the caveat that the uncertainties in the time delays 
derived from the Monte Carlo simulations described 
below may be over-estimated. 

For each round of the simulation we generate four 
sparsely-sampled fake light curves. The flux density 
for component j at a time U in the fake curves is given 
by: 

Sj(ti) = S B Rj(S (ti + Atj) + mj), j = A, B, C, D 

(11) 

where So (t) is the normalized master flux density, the 
Rj are the input flux density ratios (2.0418, 1.000, 
1.0376, and 0.3512), the Atj are the input time de- 
lays (31 d, d, 36 d, and 76 d), and is the random 
offset. The four curves are sampled with the pattern 
used in the observations (see Table 1). The flux den- 
sity error at each point in the sparsely-sampled curves 
is set to the observed flux density error for that epoch 
(see §3.4). 

The fake curves are processed in the manner de- 
scribed in §4. The best-fit time delays and flux den- 
sity ratios for each simulation are recorded. His- 
tograms of typical distributions of time delays from 
10,000 repetitions of the above procedure are shown 
in Fig. 12. The distributions are non-Gaussian, both 
in the shape of the peak of the distribution and in 
the long tail of outliers at negative delays. Thus, we 
determine the confidence limits by finding the range 
of delays inside which 95% of the simulation results 
lie, rather than fitting a Gaussian to the distribution. 
The limits are chosen such that the minimum range of 
delays that encloses both the median value and 95% 
of the simulation results is found. The results for the 
X 2 minimization and cross-correlation techniques are 
given in Table 7. We conservatively take the broader 
distribution in each case as our estimate of the 95% 
confidence contours. Thus, we estimate the time de- 
lays to be At B A = 31 ± 7 d, At B c = 36 ± 7 d, and 
AtsD = 76^10 d at 95% confidence. 

We have also run Monte Carlo simulations in which 
the sampling pattern is varied, following the method 
used in Biggs et al. (1999). The distributions of time 
delays do not differ significantly from those presented 
above. 
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Table 7 

Monte Carlo Simulation Results 











Confidence Intervals 




Quantity 


Input Value 


Method" 


68% 


90% 


95% 


AtBA 


31 


1 


28.6 - 33.0 


26.7 - 34.9 


24.4 - 37.3 


At BC 


36 


1 


33.5 - 38.7 


31.8 - 40.3 


29.0 - 43.3 


AiflD 


76 


1 


73.2 - 78.9 


71.1 - 81.0 


67.4 - 84.7 


At BA 


31 


2 


29 - 34 


27 - 36 


25 - 38 


At BO 


36 


2 


33 - 38 


31 - 40 


29 - 42 


AiflD 


76 


2 


71 - 77 


68 - 80 


66 - 82 


Sa/Sb 


2.0418 


1 


2.0386 - 2.0506 


2.0344 - 2.0549 


2.0322 - 2.0570 


Sc/Sb 


1.0376 


1 


1.0340 - 1.0401 


1.0318 - 1.0422 


1.0306 - 1.0433 


Sd / Sb 


0.3512 


1 


0.3500 - 0.3525 


0.3492 - 0.3534 


0.3487 - 0.3539 



"Method used to calculate delays or flux density ratios, 
correlation. 



minimization. 
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Fig. 10. — Composite light curve constructed as in 
Fig. 9, but using smoothed rather than raw light 
curves. The curves are smoothed with a boxcar of 
width 15 d. For clarity, only one out of every 7 points 
is shown for each curve. 



3x10 



2x10- 




Fig. 11. — Distribution of cross-correlation values ob- 
tained from 3000 Monte Carlo simulations of random- 
ized light curves. The distributions presented in this 
figure are obtained by smoothing the light curves with 
a boxcar of width 10 d. The dashed vertical lines rep- 
resent the peak cross-correlation values obtained from 
the real data. 
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6. Discussion 

The goal of monitoring a gravitational lens sys- 
tem is to measure time delays that can then be com- 
bined with a model of the lensing potential to pro- 
duce a measurement of Hq. We have been successful 
in measuring the three independent time delays in 
the B1608+656 system. Paper 2 presents a model for 
the B1608+656 system. The model is based on the 
time delays and flux density ratios presented here, and 
on positions from VLBA (Fassnacht et al. 1999) and 
HST observations of the system (Jackson et al. 1997). 
The lensing potential contains contributions from the 
two lensing objects seen in HST images of the system 
(Jackson et al. 1997), each of which is modeled as an 
elliptical isothermal mass distribution. The effects of 
varying the positions of the lensing galaxies, of chang- 
ing the nature of the lensing galaxy cores (singular 
or non-singular), and of departing from an isother- 
mal profile are all explored. The best-fit model is 
obtained through a simulated annealing process. It 
correctly reproduces the positions and flux density 
ratios of the lensed images (with the exception of the 
D/A flux density ratio). Most importantly, the pre- 
dicted time delay ratios match the observed time de- 
lay ratios to within 1%. Although the individual time 
delays depend on the Hubble constant (Ati oc /i" 1 ), 
the time delay ratios have no H dependence. Thus, 
the model for any lens system with more than two 
images must correctly reproduce the observed time 
delay ratios if it is to be used in the determination of 
Hq. In this sense, gravitational lenses which produce 
more than two images can put stronger constraints 
on lens models than can two image lenses, as long as 
the time delays can be measured. 

The B1608+656 system is the first four-image lens 
for which the three independent time delays have been 
measured and a model correctly reproduces the time 
delay ratios. The best-fit isothermal model from Pa- 
per 2 predicts time delays of AtsA = 18.0 ft." 1 d, 
At BC = 21 .4 hr x d, and At BD = 44.9 hr x d. Com- 
bining these predicted values with the observed time 
delays gives three individual determinations of Hq: 
(H )ba = 58.1 L (H )bc = 59.4, and (H )bd = 
59.1 kms -1 Mpc -1 . We combine these by calculating 
a weighted mean of (i/o)i608 = 59.0 kms~ 1 Mpc~ 1 , 
where the weights are derived from the uncertainties 
in the time delays from §5.2. Note that the B-D value 
dominates the mean since it has the smallest frac- 
tional uncertainty. In order to estimate the uncer- 



tainty on (i?o)i608j w e use the results of the Monte- 
Carlo simulation presented in §5.2. The distribution 
of time delays from the simulations are converted into 
distributions of Ho by dividing the predicted delays 
from the Paper 2 model by the simulation results. 
The resulting H distributions from the B-A, B-C, 
and B-D delays are combined into an overall distri- 
bution, which is shown in Fig. 13. The confidence 
limits are estimated by finding the ranges that en- 
close 68, 90, and 95% of the data. The resulting 95% 
confidence limits are (-ffo)i608 = 59^ kms -1 Mpc -1 . 
These confidence limits are very close to the statisti- 
cal 95% limits derived from the lens model, which in- 
cludes H as a model parameter (59lg kms -1 Mpc" 1 ; 
Paper 2). 

The above estimates of the uncertainties in the 
determination of H have not included the system- 
atic effects from the choice of the radial mass pro- 
file in the lens modeling. The estimated system- 
atic error is ±15km s _1 Mpc _1 (Paper 2). It may 
be possible to reduce this error by modeling the ex- 
tended lensed stellar emission from the background 
source, as has been done with the radio Einstein ring 
MG 1654+1346 (Kochanek 1995). This modeling ap- 
proach is being conducted by Surpi & Blandford (pri- 
vate communication). 

7. Summary 

We have presented the results of an intensive pro- 
gram of monitoring the four-image lens system B1608+656 
with the VLA. The component light curves show ^5% 
variations in flux density from which we have mea- 
sured the three independent time delays in this sys- 
tem: At B A = 31 ± 7 d, At BC = 36 ± 7 d, and 
AtBD — 76i]* d. These time delays are combined 
with the mass model of the lens presented in Paper 
2 to give H = 59tj km s -1 Mpc" 1 at 95% confi- 
dence (statistical) ±15 kms" 1 Mpc" 1 (systematic). 
The statistical uncertainties represent the 95% con- 
fidence interval. The statistical uncertainties in the 
time delays can be reduced if a stronger variation in 
the background source is observed, while the system- 
atic uncertainties may be reduced through the inclu- 
sion of the lensed extended emission in the lens mod- 
cling process. Our previous observations have shown 
that the background source in this system has varied 
by as much as 15% in the past, so we are conducting 
another program of monitoring. If a stronger vari- 
ation is detected in the new data, the uncertainties 
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Fig. 12. — Distribution of time delays recovered from 
cross-correlation analysis of 10,000 Monte Carlo sim- 
ulations of the B1608+656 component light curves. 
The distributions provide estimates of the uncertain- 
ties in the time delays. 
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Fig. 13. — Distribution of estimates of Hq from Monte 
Carlo simulations. This distribution is formed by con- 
verting the time delay simulation results into three 
H distributions, and then combining the individual 
distributions. 



on the time delays will be reduced and the accuracy 
of the measurement of H with this system will be 
improved. 
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Frazer Owen, Michael Rupen, David Rusin, Martin 
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the European Commission, TMR Program, Research 
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REFERENCES 

Bevington, P. R. 1969, Data Reduction and Error 
Analysis for the Physical Sciences, (New York: 
McGraw-Hill) 

Biggs, A. D., Browne, I. W. A., Helbig, P., Koopmans, 
L. V. E., Wilkinson, P. N., & Perley, R. A. 1999, 
MNRAS, 304, 349 

Blandford, R. & Narayan, R., 1986, 310, 568 

Blandford, R. D. & Narayan, R. 1992, ARA&A, 30, 
311 

Browne, I. W. A., Patnaik, A. R., Wilkinson, P. N., 
& Wrobel, J. M. 1998, MNRAS, 293, 257 

Fassnacht, C. D., Womble D. S., Neugebauer, G., 
Browne, I. W. A., Readhead, A. C. S., Matthews, 
K., & Pearson, T. J. 1996, ApJ, 460, L103 

Fassnacht, C. D. et al. 1999, in preparation. 

Haarsma, D. B., Hewitt, J. N., Lehar, J., & Burke, 
B. F. 1999, ApJ, 510, 64 

Hogbom, J. 1974, ApJS, 15, 417 

Impey, C. D., Falco, E. E., Kochanek, C. S., Lehar, 
J., McLeod, B. A., Rix, H.-W., Peng, C. Y., & 
Keeton, C. R. 1998, ApJ, 509, 551 



19 



Jackson, N. J., Nair, S., Browne, I. W. A. 1997, 
in Observational Cosmology with the New Radio 
Surveys, eds., M. Bremer, N. Jackson & I. Perez- 
Fournon, (Dordrecht: Kluwcr) 315 

Kochanek, C. S. 1995, ApJ, 445, 559 

Koopmans, L. V. E. & Fassnacht, C. D. 1999, ApJ, 
submitted (Paper 2) 

Kundic, T. ct al. 1995, ApJ, 455, L5 

Kundic, T. et al. 1997, ApJ, 482, 75 

Lovell, J. E. J., Jauncey, D. L., Reynolds, J. E., 
Wieringa, M. H., King, E. A., Tzioumis, A. K., 
McCulloch, P. M., & Edwards, P. G. 1998, ApJ, 
508, L51 

Moore, C. B. & Hewitt, J. N. 1997, ApJ, 491, 451 

Myers, S. T., et al. 1995, ApJ, 447, L5 

Myers, S. T., et al. 1999, in preparation 

Oscoz, A., Mediavilla, E., Goicoechea, L. J., Serra- 
Ricart, M. & Buitrago, J. 1997, ApJ, 479, L89 

Patnaik, A. R., Browne, I. W. A., Wilkinson, P. N., 
& Wrobel, J. M. 1992b, MNRAS, 254, 655 

Pelt, J., Hoff, W., Kayser, R., Refsdal, S., & 
Schramm, T. 1994, A&A, 286, 775 

Pelt, J., Kayser, R., Refsdal, S., & Schramm, T. 1996, 
A&A, 305, 97 

Pelt, J., Hjorth, J., Refsdal, S., Schild, R. & Stabell, 
R. 1998, A&A, 337, 681 

Perley, R. A. & Taylor, G. B. 1999, VLA Calibrator 
Manual 

Refsdal, S. 1964, MNRAS, 128, 307 

Schechter, P. L., et al. 1997, ApJ, 475, L85 

Schneider, P., Ehlers, J., & Falco, E. E. 1992 Gravi- 
tational Lenses, (New York: Springer- Verlag) 

Shepherd, M. C. 1997, in Astronomical Data Analysis 
Software and Systems VI, eds. G. Hunt & H. E. 
Payne, (ASP Conference Series, vl25) 77 

Snellen, I. A. G., de Bruyn, A. G., Schilizzi, R. T., 
Miley, G. K. & Myers, S. T. 1995, ApJ, 447, L9 



Walsh, D., Carswell, R. F., & Weymann, R. J. 1979, 
Nature, 279, 381 

Weymann, R. J., Latham, D., Roger, J., Angel, P., 
Green, R. F., Liebert, J. W., Turnshek, D. A., 
Turnshek, D. E., & Tyson, J. A. 1980, Nature, 285, 
641 

White, R. J. & Peterson, B. M. 1994, PASP, 106, 879 

Wilkinson, P. N., Browne, I. W. A., Patnaik, A. R., 
Wrobel, J. M., & Soratia, B. 1998, MNRAS, 300, 
790 



This 2-column preprint was prepared with the AAS IATfrjX 
macros v4.0. 



20 



